import numpy as np
import matplotlib.pyplot as plt

# 线性平移预处理，确保数据级比在可容覆盖范围
def greyModelPreprocess(dataVec):
    "Set linear-bias c for dataVec"  # 对数据进行预处理的函数

    import numpy as np  # 导入 numpy 库进行科学计算
    from scipy import io, integrate, linalg, signal  # 导入 scipy 库的相关模块
    from scipy.sparse.linalg import eigs  # 导入 sparse.linalg 模块中的 eigs 函数，用于计算特征值和特征向量
    from scipy.integrate import odeint  # 导入 integrate 模块中的 odeint 函数，用于求解常微分方程的数值解

    c = 0  # 初始化线性偏差 c
    x0 = np.array(dataVec, float)  # 将数据向量 dataVec 转化为浮点类型的 numpy 数组
    n = x0.shape[0]  # 获取数据向量的长度
    L = np.exp(-2/(n+1))  # 计算级比下限 L
    R = np.exp(2/(n+2))  # 计算级比上限 R
    xmax = x0.max()  # 获取数据向量的最大值
    xmin = x0.min()  # 获取数据向量的最小值

    if (xmin < 1):  # 如果最小值小于 1
        x0 += (1-xmin)  # 将数据向量加上一个偏差，使得最小值变为 1
        c += (1-xmin)  # 更新线性偏差 c
    xmax = x0.max()  # 更新数据向量的最大值
    xmin = x0.min()  # 更新数据向量的最小值

    lambda_ = x0[0:-1] / x0[1:]  # 计算级比，即当前元素与后一个元素的比值
    lambda_max = lambda_.max()  # 获取级比的最大值
    lambda_min = lambda_.min()  # 获取级比的最小值

    while (lambda_max > R or lambda_min < L):  # 当级比不在范围内时
        x0 += xmin  # 将数据向量加上一个偏差，使得最小值变为原最小值加上偏差
        c += xmin  # 更新线性偏差 c
        xmax = x0.max()  # 更新数据向量的最大值
        xmin = x0.min()  # 更新数据向量的最小值
        lambda_ = x0[0:-1] / x0[1:]  # 重新计算级比
        lambda_max = lambda_.max()  # 获取新的级比的最大值
        lambda_min = lambda_.min()  # 获取新的级比的最小值

    return c  # 返回线性偏差 c

# 灰色预测模型
def greyModel(dataVec, predictLen):
    "Grey Model for exponential prediction"
    import numpy as np # 导入 numpy 库进行科学计算
    from scipy import io, integrate, linalg, signal  # 导入 scipy 库的相关模块
    from scipy.sparse.linalg import eigs # 导入 sparse.linalg 模块中的 eigs 函数，用于计算特征值和特征向量
    from scipy.integrate import odeint # 导入 integrate 模块中的 odeint 函数，用于求解常微分方程的数值解

    x0 = np.array(dataVec, float) # 将数据向量 dataVec 转化为浮点类型的 numpy 数组
    n = x0.shape[0]  # 获取数据向量的长度
    x1 = np.cumsum(x0) # 计算数据向量的累加项
    B = np.array([-0.5 * (x1[0:-1] + x1[1:]), np.ones(n-1)]).T  # 构建 B 矩阵，用于计算最小二乘解析法的参数
    Y = x0[1:] # 构建 Y 向量，用于计算最小二乘解析法的参数
    u = linalg.lstsq(B, Y)[0]  # 使用最小二乘解析法求解参数 u

    def diffEqu(y, t, a, b):   # 定义差分方程，用于进行微分方程的数值求解
        return np.array(-a * y + b)

    t = np.arange(n + predictLen)   # 设置时间向量，用于微分方程的数值求解
    sol = odeint(diffEqu, x0[0], t, args=(u[0], u[1]))  # 使用 odeint 函数求解微分方程的数值解
    sol = sol.squeeze() # 压缩数组的维度，去除多余的空维度
    res = np.hstack((x0[0], np.diff(sol))) # 计算差分结果，作为预测结果
    return res   # 返回预测结果

# 输入数据
x = np.array([17730,16303, 18892, 26365, 28856, 27944, 26866, 30742, 34356])
c = greyModelPreprocess(x)
x_hat = greyModel(x+c, 2)-c

# 画图
t1 = range(x.size)
t2 = range(x_hat.size)
plt.plot(t1, x, color='r', linestyle="-", marker='*', label='True')
plt.plot(t2, x_hat, color='b', linestyle="--", marker='.', label="Predict")
plt.legend(loc='upper right')
plt.xlabel('xlabel')
plt.ylabel('ylabel')
plt.title('Prediction by Grey Model (GM(1,1))')
plt.show()

def relative_residual(dataVec, predictVec):
    "Compute the relative residual sequence"
    eps = np.abs(dataVec - predictVec) / dataVec
    return eps

def residual_test(dataVec, predictVec):
    "Perform residual test on the prediction"
    eps = relative_residual(dataVec, predictVec)
    mean_eps = np.mean(eps)
    std_eps = np.std(eps)
    return mean_eps, std_eps

mean_eps, std_eps = residual_test(x_hat, x_hat)
print("Mean relative residual:", mean_eps)
print("Standard deviation of relative residual:", std_eps)